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Motivated by recent experiments in low-dimensional trapped fermionic superfluids, we study a quasi-one- 
dimensional (quasi- ID) superfluid with a population imbalance between two hyperfine states using an exact 
mean-field solution for the order parameter When an effective "magnetic field" exceeds a critical value, the 
superfluid order parameter develops spatial inhomogeneity in the form of a soliton lattice. The soliton lattice 
generates a band of quasiparticle states inside the energy gap, which originate from the Andreev bound states lo- 
calized at the solitons. Emergence of the soliton lattice is accompanied by formation of a spin-density wave, with 
the majority fermions residing at the points in space where the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) order 
parameter vanishes. We discuss possibilities for experimental detection of the quasi- ID FFLO state using elastic 
and inelastic optical Bragg scattering and radio-frequency spectroscopy. We show that these measurements can 
provide necessary information for unambiguous identification of the spatially-inhomogeneous quasi- ID FFLO 
state and the soliton lattice formation. 

PACS numbers: 03.75.Kk, 03.75.Hh, 74.25.Gz, 74.20.Mn 



I. INTRODUCTION 

Superconducting pairing in a system with imbalanced pop- 
ulation of two fermion species (typically associated with the 
up and down spin projections) has been studied for a long 
time. Fulde and Ferrell |1], as well as Larkin and Ovchin- 
nikov [2], proposed theoretically that the superconducting or- 
der parameter in such a system is spatially inhomogeneous. 
Although several candidates for this exotic superconducting 
state (dubbed in the literature as the FFLO state) have been 
identified since then EHE), still there is no conclusive exper- 
imental evidence for the existence of this state in crystalline 
materials (see review l9l). 

Recently, great progress has been made in studying BCS- 
like pairing between fermionic neutral cold atoms. In these 
systems, it is possible to engineer a suitable attractive interac- 
tion and control population imbalance of the spin up and down 
states. Thus, cold atoms have attracted a lot of attention for 
experimental (T0lfl2) and theoretical 1 13-38] investigations 
of the FFLO state. 

The FFLO state is expected to be more stable in one di- 
mension (ID) than in three dimensions. Thus, considerable 
experimental effort was made to realize trapping potentials in 
the form of an array of ID tubes HHEL as shown in Fig.[T^. 
Although such systems have a strong ID anisotropy, it is im- 
portant to realize that they are not strictly one-dimensional be- 
cause of a non-zero tunneling amplitude t± between the tubes. 
In the quasi- ID case t± <C Ep, where Ep is the Fermi energy 
of ID motion along the tubes, the system has an open, warped 
Fermi surface shown in Fig.|T]3. Experimental results recently 
reported by the Rice group 1 12] were obtained for t± <C T, 
where T is the temperature. This regime, essentially, corre- 
sponds to an incoherent mixture of ID tubes, where each tube 
behaves independently. If, however, the strength of the inter- 
tube tunneling is increased by lowering the confining optical- 



lattice potential in the transverse directions so that t± ^> T, 
the system would cross over into the quasi- ID regime, where 
the superfluid phases on multiple tubes are locked together. 
Such a system represents a quasi- ID fermionic superfluid. 
Similar quasi- ID electronic systems have been studied exten- 
sively in solid-state physics, particularly among organic con- 
ductors and superconductors (39} [40). Quasi- ID superfluid 
states can be also realized in two and three-dimensional opti- 
cal lattices using p-orbital bands (38). 

Many theoretical papers studied strictly ID fermionic 
systems with population imbalance using analytical ap- 
proaches based on the Luttinger liquid EH [18) and the 
Bethe ansatz fT9lj23l I4TH43) . Numerical methods based on 
the density-matrix renormalization group and time-evolving 
block decimation 124^261 [28ti30l . as well as quantum Monte 
Carlo [ 31 , 32 ], were also employed. Most of these approaches 
start from the ID Luttinger liquid fixed point and take into 
account the inter- tube tunneling amplitude t± perturbatively. 
This is justified as long as t± is much smaller than the super- 
fluid transition temperature T c , see, e.g., Ref. (44). However, 
in the opposite regime t± ^> T c considered in our paper, the 
strictly ID approaches are not applicable, because the quasi- 
particle dispersion substantially deviates from the ID form. 
In this domain, the quasi- ID system is more appropriately de- 
scribed within a Fermi-liquid picture and a mean-field theory 
of superfluid pairing. 

The FFLO state in cold atoms has been recently studied in 
the quasi- ID geometry within a mean-field theory in Ref. |[T6l . 
It has been known in the solid-state literature that, within 
a mean-field Bogoliubov-de Gennes (BdG) theory, an exact 
self-consistent inhomogeneous pairing potential for a quasi- 
1D system with population imbalance has the form of a soli- 
ton lattice. The soliton-lattice solution was first obtained in 
the context of the Peierls model for charge-density waves 
in Refs. I45H48) (see reviews |49l [50)). Subsequently, this 



FIG. 1: (Color online) (a) An array of ID tubes confining cold atoms, 
where n is the tube number, and t± is the tunneling matrix element 
between adjacent tubes, (b) Fermi surfaces for the majority (f) and 
minority (D atoms in the quasi- ID limit t± <C Ef (where Ef is 
the Fermi energy of ID motion along the tube). The momentum 
components p and p± are parallel and perpendicular to the tubes, 
respectively. 



model was mapped via a particle-hole transformation onto 
the FFLO pairing problem in quasi- ID superconductors iBTT- 
l53ll . Machida and Nakanishi [52] applied these results to the 
superconducting material ErRh 4 B 4 where a strong molecu- 
lar ferromagnetic field is present. Buzdin and Polonskii l53l 
studied a similar problem for organic superconductors, where 
imbalance between the spin-up and spin-down electrons can 
be induced by an external magnetic field. An important re- 
sult found in Refs. (52] [53) is the existence of the second- 
order phase transition between the uniform and spatially- 
inhomogeneous superconducting states with an increase of the 
effective magnetic field h, which represents the difference be- 
tween the chemical potentials of the spin-up and spin-down 
electrons. When h exceeds a critical value h c , the supercon- 
ducting order parameter develops a spatially -periodic modu- 
lation in the form of a soliton lattice. 

In the present paper, we first summarize the self-consistent 
soliton-lattice solution of the BdG equations in the context of 
cold atoms and then focus on the observable physical prop- 
erties of the soliton-lattice state. In particular, we discuss 
three different spectroscopic methods for experimental de- 
tection and investigation of the soliton lattice in cold atoms. 
Optical spectroscopy of a soliton lattice for charge-density 
waves in conducting polymers was studied theoretically by 
Brazovskii and Matveenko |54|. However, because of the 
difference in coherence factors between charge-density waves 
and superconductors, there results are not directly applicable 
to the quasi- ID fermionic superfluids. Recently, there have 
been several numerical studies discussing properties of the 
FFLO state in ID geometry EUEUISllSI. However, these 
approaches are, strictly speaking, not applicable to the quasi- 
1D situation of our interest where t± ^> T c . Signatures of the 
FFLO phase have been recently studied numerically for three- 
dimensional and quasi- ID optical lattices in Refs. l35l l36ll 
which have some overlap with our results. 

The paper is organized as follows. In Sec. [TTJ we qualita- 
tively discuss a relationship between the FFLO state and the 
soliton lattice. In Sec. ITTTl we introduce the model Hamilto- 




FIG. 2: (Color online) The ID energy dispersion E(p) of the 
fermionic atoms in the hyperfine states |1) and |2) with a popula- 
tion imbalance described by the chemical potentials fi^ and /i^ and 
the Fermi momenta andp^. The third, unpopulated state |3) can 
be used for detection purposes. The atomic energy differences be- 
tween the states \2) and |1) and the states |3) and \2) are denoted as 
uj 12 and uj 23, respectively. 



nian and review basic properties of the exact solution of the 



BdG equations for the quasi- ID case. In Sec. IV we propose 
and theoretically analyze several experiments for detection of 
the soliton lattice in cold- atom settings. Finally, we conclude 
in Sec. [V] Technical details of mathematical derivations are 
relegated to Appendices. 



II. QUALITATIVE DISCUSSION 

Consider a fermionic atom, e.g., 6 Li, with two hyperfine 
states |1) and |2), e.g., F = \, mp = \ and F = \, 
mp = — \ , where F and m f are the total spin and its pro- 
jection along the quantization axis. The energies of these two 
states differ by the Zeeman splitting a; 12 due to an applied 
magnetic field, which will be not further discussed in this pa- 
per. (The Planck constant H is set to unity everywhere.) We 
denote these two states as the spin-up |1) = | t) and spin- 
down 1 2) = I 1) states. The atoms are loaded into the quasi- 
1D trap shown in Fig.[T^, where they have a parabolic disper- 
sion along the tubes and are confined in the transverse direc- 
tions. The energy dispersion of the states |1) and \2) is shown 
in Fig. [2] vs. the longitudinal momentum p along the tubes for a 
fixed transverse momentum p± . By applying radio-frequency 
(rf) radiation with the frequency uj = uj\2, it is possible to in- 
duce transitions between the states |1) and |2) and, thus, con- 
trol their relative populations. These states are long-lived, so, 
after initialization of the system, the populations of the states 
|1) and 1 2) are fixed during the time of the experiment. The 
states are characterized by the different chemical potentials p^ 
and the ID Fermi momenta p F and pp, and the densities 
of atoms per unit length p^ = p^/ir and p± = p^/ir. For 
concreteness, we assume that p^ > p^ and refer to the | t) 
and I \) states as majority and minority. The difference of the 
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FIG. 3: (Color online) Schematic plot of the quasiparticle energy 
dispersion E{p) in the inhomogeneous superfluid phase. Here £(p) 
is the normal- state dispersion, 2/x = ^ + /x^, and Q — p^ F 
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FIG. 4: (Color online) Schematic plot of the superfluid pairing for the 
LO order parameter. In each panel, the centers of the three pairs of 
arrows, representing pairs of atoms on the Fermi surface, are located 
at the point with the momentum either +Q/2 or — Q/2, which is a 
half of the Cooper pair momentum in the LO phase. 



chemical potentials is 2h — p^ — p^, where h is referred to 
as the effective magnetic field. The difference of the Fermi 
momenta Q = p F — p F results in the spin density per unit 
length n s = p^ — Pi = Q/tt. It is also useful to define the 
dimensionless spin polarization 



P = 



gt ~ Pi 



(1) 



An attractive s-wave interaction between the atoms results 
in the Cooper pairing of the fermionic atoms in the states | t) 
and ||). It is convenient to subtract the energy U12 from 
the energy of the state \2) and, thus, align the bottoms of the 
two bands, as shown in Fig. [3] In order to discuss the super- 
conducting pairing, let us make the particle-hole transforma- 
tion for the minority atoms, so that their dispersion relation 
becomes represented by the inverted parabola in Fig. [3] Be- 
cause of the mismatch of the Fermi momenta due to popula- 
tion imbalance, the conventional spatially-uniform BCS pair- 
ing potential is not favorable, since it would open an energy 
gap away from the chemical potentials of the atoms. Larkin 
and Ovchinikov (LO) |2] proposed that the pairing poten- 
tial should be non-uniform and have the spatial dependence 
A(x) oc s'm(Qx) l55l . This order parameter couples fermions 
having the difference ±Q of the momenta \p^\ and \p±\, so it 
opens energy gaps at the chemical potentials for both major- 
ity and minority atoms, as shown in Fig. [3] Then, the lower 
band is populated by the Bogoliubov quasiparticles with both 
spins up and down, the middle band is populated only by the 
quasiparticles with spin up, and the upper band is empty. 

Fig. [3] only illustrates the ID dependence of the quasiparti- 
cle energy E(jp) on the longitudinal momentum p in the LO 
state. However, even for a non-zero t±, the LO order parame- 
ter A(x) oc sin(Qx) still opens a gap on the whole quasi- ID 
Fermi surface, as illustrated in Fig. [4] and explained below. 
For simplicity, let us consider the tunneling amplitude t± be- 
tween the tubes only in one transverse direction, as shown in 
Fig.[T^, and denote the corresponding momentum p±. Then, 
the normal-state quasi- ID energy dispersion is 



e*(P:P±.) 



2 in 



2t± cos(p_\_d) — p — crh, (2) 



where m is the mass of an atom, d is the inter-tube spacing, 
and ah = ±h for a =t,l- In the case where t± : h <C p, 
we can linearize the energy dispersion relation along the tube 
direction and obtain 



£*(P,P±) ~ ±v f (ptPf) -2t± cos(p±d), 



(3) 



where v F = pp/m is the Fermi velocity, p a F = Pf + crQ/2, 

and Q = 2h/vp. The Fermi surfaces e a (p,p_\_) = obtained 

from Eq. ^ are shown in Fig.[TJ). 

As illustrated in Fig. HI the majority and minority fermions 
■f i — 

on the +p l F and — p F branches of the Fermi surface pair in 
such a manner that the total momentum of each pair is +Q 
and does not dependent of the transverse momenta ±p±_ of 
the paired fermions. Correspondingly, the total momentum 
is — Q for each pair on the — p F and +p F branches. So, 
the LO order parameter with the two momentum components 
A(x) oc (e l ® x — e~ tQx ) opens an energy gap on the whole 
Fermi surface for all values of p±. Thus, one expects to have 
a stable LO phase in the quasi- ID geometry. However, this 
result is valid only for the linearized energy dispersion ^ 
and does not apply for a generic three-dimensional dispersion, 
e.g., for a spherical Fermi surface. 

The linearized dispersion ([3| can be rewritten in the form 



±v f (pTPf) ~ crh, 



where 



Pf(pa 



- p F 



2t± 

v F 



cos(p±d) 



(4) 



(5) 



is the longitudinal Fermi momentum as a function the trans- 
verse momentum p± at h = 0. As will be shown in Sec.llllj it 



is possible to make a gauge transformation of the fermion op- 
erators and eliminate pp along with any explicit dependence 
on p± and t± from the Hamiltonian of the problem. Physi- 
cally, this gauge transformation corresponds to measuring the 
longitudinal momentum p from the local Fermi momentum 
Pf(p±) at each point on the quasi- ID Fermi surface. After 
this transformation, the mathematical problem formally be- 
comes one-dimensional, and an exact solution of the mean- 
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FIG. 5: (Color online) Schematic plots of the density of states (top 
row), superconducting order parameter A(x) (middle row), and the 
spin density p s {x) (bottom row) for different values of the effective 
magnetic field h. (a) For h < h c , the ground state of the system 
is uniform, as in the BCS theory, (b) For one unpaired atom h = 
h c + Sh, A(x) has a kink soliton. (c) For a small spin imbalance 
h > h c , A(x) forms a soliton lattice, (d) For a large imbalance 
h ^> h c corresponding to the LO limit, A(x) is sinusoidal. Notice 
that p s (x) peaks at the points in space where A(x) vanishes. 



field equations can be obtained, which reduces to the LO or- 
der parameter and the soliton lattice in different limits. Nev- 
ertheless, it is important to remember that, although t± can 
be formally eliminated from the mean-field Hamiltonian, the 
physical problem remains quasi- ID, and the presence of t± 
stabilizes fluctuations of the order parameter and makes the 
mean-field approach applicable. 

Fig. [3] illustrates the case where h ^> A. When A and 
h are comparable, the two energy gaps in Fig. [3] are coupled 
and cannot be treated independently. The problem becomes 
mathematically complicated in this case. Fortunately, an exact 
self-consistent solution of the BdG equation in the quasi- ID 
case can be obtained (45l[52j[53l. The exact solution shows 
that, in general, A(x) is given by a periodic Jacobi elliptic 
function, such that it reduces to sin(Qx) in the limit h ^> A. 
The exact A(x) represents the so-called finite-zone potential, 
which opens only two gaps in Fig.|3]at the chemical potentials 
/i^ and nowhere else. The exact solution shows that there is 
a critical value h c = 2A /7r such that the order parameter is 
uniform A(x) = const = Ao for h < h c , where Ao is the 
value of the BCS superconducting gap for h = 0. 

For h slightly above h c , the order parameter develops a se- 
ries of kink solitons, where A(x) changes sign across each 
kink soliton. In order to explain soliton formation qualita- 
tively, let us consider the sequence of states shown in Fig. [5] 
In this figure, the middle row shows the spatial dependence of 
A(x), the top row shows the density of quasiparticle states, 
and the bottom row shows the spatial dependence of the lo- 
cal spin density p s {x). Panel (a) shows that, for h < h c , the 



system has the conventional BCS order parameter with the 
uniform A(x) = Ao, a single gap in the density of states, 
and zero spin density p s = 0, i.e., no spin imbalance. Now, 
suppose we add one extra majority atom to the system. Be- 
cause this atom cannot form a pair, it would have to pop- 
ulate an energy level above the gap, which would cost the 
energy A . However, if the system develops a kink, i.e., a 
sign change of A(x), as shown in Panel (b), this configuration 
has an Andreev bound state in the middle of the gap, which 
can be occupied by the extra atom. The total energy cost for 
creation of the soliton and occupation of the midgap state is 
(2/7r)Ao, which is lower than Ao for the uniform configura- 
tion. Thus, the system spontaneous creates the kink soliton 
shown in Panel (b). The spin density p s (x) of the extra atom 
is concentrated near the soliton. Soliton creation has been dis- 
cussed in the literature for charge-density waves |45j|46l, for 
mesoscopic superconducting wires |[56lL and in the cold-atom 
context l57l . 

When a few atoms with spins up are added to the system, 
each atom creates a kink soliton. The solitons repel each other 
and form a periodically spaced soliton lattice shown in Panel 
(c). The order parameter A(x) experiences a series of sign 
changes. The midgap states from different solitons hybridize 
and form a band in the middle of the main energy gap. The 
spin density p s {x) consists of a series of spikes originating 
from each soliton. When many atoms with spins up are added, 
A(x) becomes sinusoidal, corresponding to the LO limit, as 
shown in Panel (d). The midgap band expands and occupies 
most of the former energy gap, leaving two small gaps above 
and below, in agreement with Fig. [5] The spin density p s (x) is 
sinusoidally modulated with a small amplitude and the wave- 
length a half of that for A(x). 

In the next Section, we present a detailed mathematical 
derivation of the results qualitatively discussed above. Then, 



in Sec. IV we discuss three spectroscopic techniques for ex- 
perimental detection of the soliton lattice. 



III. THEORETICAL MODEL AND EXACT SOLUTION 

Let us consider a two-dimensional (2D) array of parallel ID 
tubes in the x direction with the inter-tube spacings d y and d z 
in the y and z directions, as sketched in Fig.[T^. In the second- 
quantized formalism, the single-particle Hamiltonian for this 
quasi- ID system is 



Ho = ^2 dx ^tiA x ) 

nr. t J 



d 2 
'2m 



/ dxt± 

( 71,71/) ,<t 



$L*( X )$n'A X ) +H.C. 



, (6) 



where ipn,a(%) is the fermion annihilation operator for an 
atom in the state a =t, I on the tube n, where n = (n y , n z ) 
is the 2D index of the tube, and t± the amplitude of tunneling 
between adjacent tubes. An attractive interaction between the 
atoms leads to the BCS pairing, which, at the mean-field level, 
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is described by the Hamiltonian 



Hsc : 



-E 



dx 



t 



(7) 



where the pairing potential is determined self-consistently 

A n (x) = g (^ n ^(x)^ n ^(x)) (8) 

with the s-wave interaction amplitude g. We assume that pair- 
ing is local in the real space, and the inter-tube pairing po- 
tential can be neglected. Given that we consider identical 
tubes, we take the pairing potential to be independent of n, 
i.e., A n (x) = A(x). We study the problem at zero tempera- 
ture T = 0, so the brackets in Eq. ([5]) represent averaging with 
respect to the ground state. 

After the Fourier transform Vv(^P±) = 
X^n, e ~ ip± ' p ' n " l Pn,a(z)> where p± = (p y ,Pz) is the transverse 
momentum, and p n = (n y d yj n z d z ) is the 2D vector in the 
(y,z) plane, the full Hamiltonian H = Hq + Hsc becomes 



H = ^2 / dx^l(x,p±) 



^+C(p J .)-M-^ 



y^ldx A(x)^(x,p±)^(x,—p±) + B..c. 



4>*(x,p±) 

(9) 



Here^(p^) = —2t±[cos(p y d y )+cos(p z d z )] is the transverse 
dispersion, and the sum over the transverse momenta means 
E PX = d y d z J^ /dy dp y J^ /dz dp z /(27r) 2 . The Hamilto- 
nian |9j) is the sum H = ^Z p± #id(p_l) of ID Hamiltonians 
Hid with fixed values of p± 

/r q2 
dx$l(x) 



2m 

dx \A(x) ?pt(x) $\(x) 



ah 



H.c. 



(10) 



where jl(p±) = [i — is the renormalized chemical po- 

tential, and the argument p± in ^(x) is implied. The self- 
consistency condition (ISl) has the form 



A 0) =9 z2(^i(x,p±)^(x,-p±)). 



p± 



(id 



The integrand in the Hamiltonian ( fT0| ) can be written in the 
matrix form as 



$\(x) 



2m 



h 



A(x) 



2fjj(x) 

^\(x) 



(12) 



where the overline in A denotes complex conjugation. The 
Hamiltonian Hid also acquires a constant term ^2 p (p 2 /2m — 
jl + h) originating from the fermion commutation relation 
ipfyl = 1 — The matrix |l2| ) can be diagonalized by 

the Bogoliubov transformation 

X r- n (I 3 ) 



where the sums are taken over the eigenstates A specified be- 
low. The coefficients U\(x) and are selected so that the 
Nambu spinor &\(x) — [U\(x), V\(x)] T is an eigenstate of 
the BdG equation with an eigenvalue E\ 



(14) 



We also impose the normalization condition 

J L L dx[\U x (x)\ 2 + \V x (x)\ 2 ] = 1, (15) 



where 2L is the length of the system. The BdG equation ( 14 ) 
has the particle-hole symmetry: if [U\(x), V\(x)] is an eigen- 
state with the energy E\, then [— V\(x), U\(x),] is also an 
eigenstate with the energy — Ey. 



U x (x) 
V x (x) 



-V x (xj 
U x (x) _ 



E x 



-E x . 



(16) 



This property allows us to restrict summations over eigenen- 
ergies in subsequent calculations only to the positive values 
£a>0. 

Substituting Eq. ( [T3] ) into Eq. ( [T2| ) and utilizing the prop- 
erties of U\(x) and V\(x) outlined above, we diagonalize the 
Hamiltonian Hid 

Hid = Yl E - oh) t^Ta,. + £b- (17) 

S A >0(7=U 



Here £0 is the reference energy given by 

£0 = - £ £a + ££ 



(18) 



where the first sum originates from the fermion commutation 



relation 7a,+7a 



1 



■ 7a,;7a,; 

o2 



and 



P 
2m 



(19) 



is the normal-state dispersion relation. The occupation num- 
bers of the single-particle states for the Hamiltonian ( p7| ) are 



(7L7A,a) =n F (E x - ah), 



(20) 



where np{E) is the Fermi distribution function, which re- 
duces to the step function n F (E) = 6{-E) at T = 0. The 
ground state energy of the system is I52ll56l 



Li 

Fo=£o+ Yl (E\ - vh)n F (E x - ah) + j 

tp \ r\ _ 



A 2 (x) 



(21) 



where the last term originates from the Hubbard-Stratonovich 
transformation of the interaction between the fermions. 
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The pairing potential A(x) in Eq. fl4| ) can be selected to be 
real. Now we linearize the dispersion relation in the longitudi- 
nal direction near the Fermi surface and use the approximation 
similar to Eq. ([4} in the BdG equation 



v F [i(-l) a d x -p F ] 
A(x) 

E x *{ a \x), 



A Or) 
-v F [i(-l) a d x -Pf 



where 



p F = ^2m]± w p F - £(p±)/v F 



(22) 



(23) 



is the renormalized Fermi momentum, and the index a — 
1,2 corresponds to the right- and left-moving atoms. We 
seek solutions of the BdG equation in the form ^\(x) — 

E** { x\x) with 



^W-f^^l^f^if) * %PFX 
A v ; \Vx,i{x) J \v x {x) 



(<*=2)/ x = fUxA X )\ _ ( V \( X ) 



(24) 



In Eq. ( [24] ), we used the symmetry between the a = 1 and 
2 components: u\ = ^a,i = ^a,2 and = ^a,i = ^a,2- 
Substituting Eq. ( [24] ) into Eq. ( [22] ), we eliminate j5p (p^ ) from 
the BdG equation and obtain the ID equation for the slowly 
varying envelope functions u x (x) and v\ (x) 



-iv F d x A(x) 
A(x) iv F d x 



u x(x)\ _ 
v\(x) 



u\(x) 
v\(x) 



(25) 



In terms of the amplitudes u\ and v\, the self-consistency 
condition ( [TT] ) now reads 

A(x) = g ^2[u x (x)v x (x) + u\(x)v x (x)} 

a (26) 
x [l-n F (E x + h)-n F (E x -h)]. 

Summation over p± in Eq. ( [26] ) is omitted, because solutions 
25] ) do not depend on p± , and the phase factors from 
containing p F (p±) cancel out in Eq. (26). Thus, as a 



result of the linearization approximation, we managed to elim- 
inate the transverse attributes t± and p± from the BdG equa- 
tion ( [25] ) and the self-consistency condition ( [26] ) and reduce 
the ground-state mean-field problem to a purely ID formula- 
tion, effectively corresponding to a single tube. Deviations of 
the actual dispersion relation from the linearization approxi- 
mation are of the order of t\/E F , so the approximation em- 
ployed in this paper is applicable when t 2 ± /E F <C T c <C 
t ± <^E F . 

Although t± has been eliminated from calculation of the 
mean-field ground state, a non-zero value of t± is crucially 
important when considering fluctuations of the order param- 
eter near the ground-state configuration. While the ground 
state A(x) depends only on the coordinate x, the fluctuating 
order parameter depends on all three coordinates A(x,y,z). 



A non-zero tunneling amplitude t± produces transverse phase 
stiffness in the effective action, which is proportional to 
t 2 ± [(d y <f>) 2 + (d z &) 2 } with <3> being the phase of the super- 
conducting order parameter. As a result, phase fluctuations 
are suppressed due to the three-dimensional anisotropic stiff- 
ness (as opposed to ID stiffness for uncoupled chains), so the 
true long-range order is stabilized at low enough temperatures, 
and the mean-field approximation is justified. Stability of the 
FFLO phase in quasi- ID geometry was recently confirmed 
numerically in Ref. l36l . 

A self-consistent solution of Eqs. ( [25] ) and ( [26] ) was first de- 
rived in Refs. l33[46l in the context of charge-density waves 
and than subsequently extended to various other systems J47J 
l48l . including inhomogeneous superconductors 15TH531 . To 
keep our discussion self-contained, we outline the properties 
of the solution and refer the reader to the literature for further 
details. Let us introduce the functions 



f ± -- 



;(u\±iV\), 

and rewrite Eq. d25|» in the following form: 



A 
it o 



ft 
fx 



ft 

fx 



(27) 



(28) 



where A = i[—vpd x + A (a;)] and A< = —i[vpd x + A (a;)] 
are mutually adjoint operators. After simple manipulations, 
Eq. (|28|> can be written in the supersymmetric (SUSY) form 



(Aft 
1 iU 



ft 

fx 



Ei 



ft 
fx 



(29) 



Equation ( |29l > corresponds to the N = 2 SUSY quantum me- 
chanics introduced by Witten l58l . where A(x) is the SUSY 
potential. This connection will be important for a discussion 
of zero-energy bound states at domain walls. In the explicit 



form, Eq. (129b reads 

, d 2 



dx 2 



Ei-A z (x)±v F 



dA(x) 
dx 



f£(x) = 0, (30) 



while the self-consistency equation ( 26 ) is 

= g^^r [v F d x + 2A(x)] \f+(x)\ 2 

x [l-n F (E x + h)-n F (E x -h)]. 



(31) 



The Schrodinger equations ( [30] ), with the effective potentials 
determined by the order parameter A(x), are integrable for a 
special type of reflectionless potentials. A solution for A(x) 
minimizing the ground-state energy To can be written in terms 
of the Jacobi elliptic functions specified by the modulus k 



A(s) 
A 2 

where 



= (l-fc')sn[(l + fc')C,H 



, 2 sn(C, fc)cn(C, k) 
' dn(C,fc) : 



C = x- 



k! = x/l - k 2 , 



k! 



(32) 



(33) 
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The second equality in Eq. ( [32] ) follows from the properties of 
the elliptic functions. The parameter A 2 and the modulus k 
of the Jacobi elliptic function have to be determined from the 
minimization of the ground-state energy in Eq. ( [2T] ), see 
Ref. l52l for more details. The period Z of the order parameter 
A(x) is 



(34) 



where K(k) is the complete elliptic integral of the first kind. 
Upon substitution of Eq. ([32]) into the ground-state energy 
and minimization with respect to A 2 and k, one finds B31I521 



A 2 7T 



(35) 



where Ao is the BCS energy gap for a homogeneous unpo- 
larized system at h = 0, and E(fc) is the complete elliptic 
integral of the second kind. The parameter Ao is a convenient 
way to characterize the strength of the attractive interaction g 
in the system. Eliminating A 2 from Eq. ( [35] ), we obtain an 
equation for the modulus k in terms of the given values of the 
effective magnetic field h and the BCS gap Aq: 



k 



E(k) 



2 Ao 

7T h 



(36) 



The values of k are restricted to k < 1 because of the proper- 
ties of the elliptic functions. With the increase of h, Eq. ([36]) 
acquires a solution at the critical value h c = 2Ao/7r, where 
k —> 1 and E(fc = 1) = 1. At h = h c , the solution A(x) cor- 
responds to a single soliton, as shown in Panel (b) in Fig. [5] 
With the further increase of h, the value of k given by Eq. ( |36| ) 
decreases, so the soliton lattice period Z in Eq. ( [34] ) also de- 
creases, corresponding to Panel (c) in Fig. [5] In the LO limit 
h ^ Aq, Eq. ( [36] ) gives k <C 1, and the period Z in Eq. ( [36] ) 
becomes Z = nvp/h. In this limit, the order parameter takes 
the LO form A(x) « (A 2 /c 2 /2) sin(27nr/Z). 

The solution for the amplitudes f ± (x) can be expressed in 
terms of the function 7(C) satisfying the following equation 



A 2 



= 4 7 ( 7 -f 2 2 )(f|- 7 ), 



(37) 



where the variable ( is defined in Eq. ( [33] ). Here we intro- 
duced the parameters £3 = A 2 and £ 2 = A 2 /V. The parame- 
ters £ 2 and £3 are the band edges in the single-particle excita- 
tion spectrum shown in Fig. [6] The spectrum in Fig. [^corre- 
sponds to the spatially-inhomogenous order parameter A(x) 
in Eq. ( [32] ) and represents a zoom-in of the spectrum shown in 
Fig. [3] in the vicinity of the Fermi momentum p « pp. 
The solution of Eq. ( [37] ) has the form 



2 7 (C)=^+^ 2 -A 2 (C) + A 2 



(38) 



One can check by direct substitution that the solution of 



Eq. ( 30 ) reads 



fU x ) 



2LA> 



exp 



ib 
vf 



fRldx' 



70') 



(39) 



where R X = E 2 X (E 2 X 



2 )y±J X — f |), 2L is the length of the 



system, and A\ is a normalization factor 



^4 a = E 



2 f2 E(fe) 



(40) 



The index b = ± in Eq. ( [39] ) distinguishes quasiparticles re- 
siding on the two energy branches of the spectrum shown in 
Fig. [6] and allows us to define unambiguously the momentum 
Pa, b for a given value of E\. Using Eqs. ( [25] ) and ( [27] ), we 
find the other solution 

fx, b ( x ) = cx, b e-^f+ b (x + l/2), (41) 
where the coefficient and the phase ip\^ are 



c x , b = sign(bE x ) x 

^A,6 



b 



+1 

R x dy 



1, |£a|<£ 2 , 



(42) 







7(2/)' 



We emphasize that, in order to recover correct amplitudes 
and v\^, it is necessary to know the relative phase be- 
tween fxb( x ) anc * f\b( x ) m Eq. This su bject was not 
discussed in the previous work focusing on thermodynam- 
ics of the soliton-lattice state (45] [47] [5TV[53lL because the 
order parameter A(x) ( [3T] ) depends only on and 
the ground- state energy can be calculated without invoking 
the relative phase. However, in general, response functions 
involve matrix elements between various Bogoliubov ampli- 
tudes, and the relative phase is absolutely necessary for main- 
taining the particle-hole symmetry. Note that the solution 



given by Eqs. ( [32] ) and ( [39] ) is a particular case of a more gen- 
eral type of multi-periodic solutions |59l . 

Because the pairing potential A(x) is periodic, the am- 



plitudes fxb( x ) m Eq. (|39| must satisfy the Bloch theorem 



a iP\,bl f± 



f± b {x + I) = e 
Here the subscript b 
Pa, 6 corresponding to a given eigenenergy E\. Thus, we can 
write 



Xh (x) with the quasi-momenta p\^. 
± labels the two momentum branches 



Jp\, b XA± 



^a,6<»> 



(43) 



where (j) x b [x + Z) = 0^ 5 (x) is a periodic function, which can 
be expanded in the Fourier series 

00 

<t>U*) = E ^He"" 91 , 

m— — 00 



2tt 

T" 



(44) 



Using the energy-momentum relation p\^ = Pb{E\\ we can 
obtain the energy dispersion E\(p\) shown in Fig. [6] and the 
density of states per spin Pb(E\) — \dp\^/ dE\\/2ir. Below 
we omit the eigenstate label A for brevity. First we consider 
the energy range < E < £ 2 and find from Eq. ([39]) H71I521 : 



Pb(E) = 



bA 2 
~v F K(k) 



M(<p E ,k'), E<£ 2 



(45) 
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FIG. 6: (Color online) The single-particle energy dispersion relation 
E(p) in the quasi- ID FFLO state described by Eq. ( [32] ). The spec- 
trum is obtained by solving Eqs. ( |45] > and < [47| ) for given values of 
the momentum p and the parameter b = =L This plot is a zoom into 
Fig.[3]near the Fermi momentum p^. 



where cfE = arcsin(|J5|/^2) and 

.M(<^, fc 7 ) = [E(fc) - K(fc)]F(^, fc') + K(fc)£(y>tf, /c 7 ) 

(46) 

with F((pE,k f ) and E((p F ,k f ) being the incomplete elliptic 
integrals of the first and the second kinds, correspondingly. 
The function A4((pE,k') remains positive for all values of the 
energy E. The momenta corresponding to the band edge en- 
ergy £ 2 = A 2 k' arep(£ 2 ) = ±(7rA 2 /2v F )K(k) = ±Q/2, as 
shown in Fig. [6] 

Similarly, for the energies above the upper band edge E > 
S 3 , we get 



Pb(E) 



b_ 

v F 



(E*-Si)(E* 



K(fc) 



E 2 



(47) 



E>6^ 



where (p F = aicsin(Es/\E\). Notice the minus sign in front 
of the square brackets in Eq. ( [47] ), which indicates that the 
branch index b actually changes across the gap, as shown 
in Fig. [6] For the energies away from the gap E <C £ 2 
and E ^> £3, Eqs. ( [45) and ( [47] ) reproduce the original lin- 
earized dispersion relation E = v F p in Eq. ([4). [The quasi- 
momentum p, introduced in Eq. ( [43] ) and used in the rest of 
the paper, actually corresponds to p — p F , where the latter 
momentum p is the one used in Sec. [II]] Equations ( [45] ) and 
( [47] ) determine the single-particle spectrum of the quasi- ID 
FFLO state, which is shown in Fig. [6] 

The single-particle density of states p(E) = \dp/dE\/27r 
is given by the following expression I 47l 



p(E) = p F 



\2E 2 



£l 



£1 



(A 2 



V(£ 2 -£f)(£ 2 -£!) 



PF 



2ttvf ' 



(48) 



where p F is the density of states in the normal state and 



1 

(A 2 ) = \ J A 2 (z) dx = A 2 ( 2 - fc 2 - 2 |||) . (49) 



Plots of the density of states for difference values of h are 
shown in the top row in Fig. [5] At the point of transition into 
the spatially homogeneous superfluid phase (k —> 1 and h — >> 
h c ), Eq. (48 ) reproduces the BCS density of states. 

Having obtained the exact Bogoliubov amplitudes, we can 
compute the spin density p s (x): 



p s (x) = (s z (x)) = {^\(x)^(x)) - (^l(x)^(x)), (50) 

where s z (x) is the spin-density operator. Using Eqs. fl3] ) and 
([20]), we find 



(1>l(x)il> a (x)) = ^{\U x (x)\ 2 n F (E x - ah) 

A 

+ \V x {x)\ 2 [l-n F {E x + ah)]} 



(51) 



and 



p s (x) = ^[\U x (x)\ 2 + \V x (x)\ 2 ]x 

a (52) 
x [n F (E x -h)-n F (E x + h)]. 

where the sum is taken over positive E\, as discussed below. 
Eq. ( [52] ) has a clear physical interpretation. First, at T = 0, the 
contribution to the spin density comes only from the middle 
band with the energies \E X \ < h, which is occupied by the 
majority spin only, as shown in Fig. [6] Second, for a given 
value of A, the spatial integral of the spin density is equal 
to one because of the normalization condition for the Nambu 
spinors ( [T5] ). Thus, at T = 0, the integral 

Us = jl I dx p*w = ^zJ2°( h - E ^ (53) 

J-L x 

/oo /< 
dEp(E)9(h-\E\) = 4 / 
-00 Jo 



<^' 2 dp = Q 

2lT 7T 



relates the spin imbalance n s and wavevector Q = p F — p F , 
as discussed in Sec. [II] The factor of 2 in the second line of 
Eq. ( [53] ) comes from summation over the index a representing 
the ±pf branches, see Fig. [3] 

It is instructive to consider the limit h —> h c , where the 
imbalance corresponds to just one excessive majority atom 
n s = 1/2L. In this regime, the order parameter forms a do- 
main wall A(x) = A tanh(A /^F^), which binds the un- 
paired atom. There is a single normalizable midgap state at 
E x = 0, which is localized at the domain wall and is charac- 
terized by the Witten topological index l58l 



U (x) 
Vo(x) 



1 



2v f cosh 



O) 



cos(p F x) 
i sm(p F x) 



(54) 
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Using Eq. ( [54] ), we obtain the spin density from Eq. ([52]) 
A 1 



Ps(x) 



2v F cosh 2 (xA /^F)' 



(55) 



which, indeed, corresponds to the unpaired spin localized at 
the domain wall around x = 0. With the increase of the pop- 
ulation imbalance, it becomes energetically favorable for the 
system to create more domain walls in order to accommodate 
excessive spins. This regime corresponds to the soliton lattice. 
In general, we can calculate p s (x) for an arbitrary population 
imbalance using the exact solution discussed above. Using 
Eqs. ([24]) and p7| ), we rewrite Eq. ([52]) at T = as 



psix) = E(I/a + (*)I 2 + I/a~(*)I 2 ) 0(h - E X ), 

A 



(56) 



where we omitted the fast-oscillating terms with the wavevec- 
tor 2pp. After substituting Eqs. ( [39] ) and ( |4T] ) into Eq. ( [56] ) 
and using the identity ^ A f(E x ) = 2L j dE p(E) f(E), we 
find 



Psjx) _ 

= E(fc') 



r^2 
JO 



dE 2E 2 



f 2 



■A 2 (x) 

2A 2 ^/^ 2 -E 2 \£ 2 -E 2 ) 



(57) 



l-k r 



■K(k f ){l + (l-k^)sn 2 [(l+k')CM}, 



The plots of A(x) and p s (x) are shown in Fig. [7] In the 
soliton-lattice limit depicted in Panel (a), the spin-density 
spikes are well pronounced, because the distance between 
solitons is fairly long. With the increase of the spin imbalance, 
the solitons start to overlap, and the system crosses over to the 
LO phase with the sinusoidal A(x) shown in Panel (b). The 
midgap states become hybridized and extended over many do- 
main walls. As a result, the amplitude of spin-density modu- 
lation becomes small compared with the uniform spin back- 
ground, as shown in Panel (b). Thus, it may be difficult to 
detect the modulation of p s (x) in the LO limit experimentally. 

Finally, we point out that the spin density p 3 (x) = p s (x + 
1/2) is a periodic function with the period 1/2 = ir/Q, so 
it can be expanded into a Fourier series with the coefficients 
p s {m) 



p s (x) = 2 Ps( m ) cos(2Qmx), 



p s (m 



m=0 

4 r l/4 



dx p s (x) cos(2Qmx). 



(58) 



(59) 



Because of the symmetry p s (x) — p s (—x), the Fourier expan- 
sion has only cosine functions. The coefficients p s (m) can be 
computed using Eq. ([57 ) 

E(fc') - ^-^K(fc')l 5^o 



Ps(m) 



4p F A 2 



(1 - k' 2 ) 2 



■K(k')Y(k',m), 



(60) 



where the dimensionless function Y(k', m) is defined as 

Y(k',m)= / dz sn 2 \K(y)z, v\ cos(7rmz) (61) 
Jo 
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FIG. 7: (Color online) The pairing potential A(x) and the spin den- 
sity p s (x) calculated for the quasi- ID FFLO state, (a) The soliton- 
lattice regime for h close to the critical value h c . The spikes of p s (x) 
are aligned with the kink solitons, where A(x) changes sign, (b) The 
Larkin-Ovchinnikov regime for h ^> A . The spin-density modula- 
tion is small compared with the uniform spin background. 



with v = {1 — k f )/(l+k f ). The Fourier coefficients p s (m) can 
be calculated by expressing the sn function in terms of an infi- 
nite series and taking the spatial integral over z, see Ref. l52l . 
We will use the expansion ( [58] ) in the next Section when dis- 
cussing the elastic optical Bragg scattering experiments. 

Here we conclude our overview of the exact mean-field 
pairing solution for a quasi- ID Fermi system with a spin- 
population imbalance. One of the hallmarks of this solution is 
the midgap energy band populated by the majority spins. In 
the next Section, we propose and theoretically analyze several 
experiments for detection of the soliton lattice. 



IV. EXPERIMENTAL DETECTION OF THE SOLITON 
LATTICE 

The soliton lattice in the quasi- ID FFLO state can be 
detected by various experimental techniques, such as po- 
larization phase-contrast imaging (12] |60), Bragg diffrac- 
tion (6HE2), radio-frequency spectroscopy |63, 64], quantum 
polarization spectroscopy l65l and quantum spin-noise spec- 
troscopy (66] [67). Different experimental techniques have 
advantages and disadvantages. In particular, the polariza- 
tion phase-contrast imaging used in the recent Rice experi- 
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ment 1 12] does not have sufficient spatial resolution to resolve 
modulation of p s {x). Thus, alternative experimental probes 
are needed to identify the FFLO state unambiguously. Also, 
it is important to map out the parameter space where the sig- 
natures of the FFLO state are most prominent. In this Sec- 
tion, we propose and analyze theoretically three experimental 
approaches which may be the most suitable for a clear de- 
tection of the FFLO phase. Specifically, we consider the opti- 
cal Bragg diffraction on spin-density modulation, the inelastic 
Bragg scattering, and the radio-frequency spectroscopy. 

A. Optical Bragg diffraction on spin-density modulation 

Spatial modulations of the spin density p s (x) in the FFLO 
state can be detected using the spin-dependent Bragg scatter- 
ing of light. This method is a cold-atom analog of the elastic 
polarized-neutron scattering widely used in condensed-matter 
physics, particularly for detecting magnetism and determin- 
ing the symmetry of a superconducting order parameter, see, 
e.g., Ref. [68]. The optical elastic Bragg scattering was pro- 
posed in Ref. l62l for detection of antiferromagnetism in cold 
atoms described by the fermionic Hubbard model, but it can 
be also adapted for our problem. This method is based on 
the following observation. If the frequency uj of incident 
light is tuned halfway between the energy distance from the 
level 1 3) to the two hyperfine levels |1) and |2) in Fig. [5J 
i.e., uj = U23 + CJ12/2, then the light couples to the popu- 
lation imbalance of the levels |1) and |2), i.e., to the local 
spin density p s (x) = p^{x) — Pi(x). As explained in more 
detail in Appendix [A| this happens because the scattering ma- 
trix elements of light on the atoms in the hyperfine states |1) 
and 1 2) have opposite signs, so they produce opposite phase 
shifts for the light. Because the light frequency is detuned 
from atomic transitions in this regime, there is no photon ab- 
sorption, whereas photon scattering is maximally sensitive to 
the spin-density modulation. Thus, the cross section of elas- 
tic scattering is proportional to the static spin- structure fac- 
tor S(q x ) 1 62], which is determined by the Fourier transform 
Ps(q x ) of the spin density p 3 (x): S(q x ) oc \p s (q x )\ 2 . Be- 
low, we derive this relation explicitly for the three-band model 
shown in Fig. [2] 

The energy dispersion of an atom in the state |3) is es(p) = 
E3 -\-p 2 /2m, where E3 = 0023 + 0J12 is counted from the bot- 
tom of the lowest band in Fig. [2] We assume that the state |3) 
is not populated (i.e., ps — 0), and the interaction amplitudes 
between atoms in the state |3) and with atoms in the other 
states |1) and |2) are negligible. Let us introduce the annihi- 
lation operators ^s(x) for an atom in the state |3) and a k for 
a photon with momentum k. The interaction of light with the 
atoms is governed by the Hamiltonian H int = (H a + H\ ) , 
where the operator H a describes atomic transitions from the 
state a to the state |3) with absorption of a photon 

r L 

H a = Y, dxT k e- ik ^a k ^l(x)^ a (x). (62) 

h J —L 

The amplitudes T k contain microscopic information about 




FIG. 8: (Color online) Schematic plot of the elastic Bragg diffraction 
experiment. Photons with the momentum \ k\ — 2tt/\ are elastically 
scattered to the momentum k' by the spin-density modulation with 
the period 1/2. The Bragg condition for constructive interference is 
I sin = mX (or 2k x = 2Qm), where m is an integer. 

atomic transitions, such as dipole matrix elements and light 
polarization. To simplify notation, Eq. ( [62] ) is written for one 
tube, but summation over all tubes is implicitly assumed. 

We consider a process where an incoming photon with the 
momentum k = (k x ,k±) and energy uj k scatters into the 
state with the momentum k! — (—k x ,k±) and energy uj k ', 
as shown in Fig. [8] Since \k\ = \k f \, the photon energies 
= Wk> are the same, so the scattering process is elastic. At 
the same time, the momentum transfer q x = k x — k' x = 2k x 
from the photon to the atoms in the x direction allows one to 
probe the spatial spin-density modulation in the FFLO state. 

The transition amplitude M k , k for this elastic process is 

given by the second-order perturbation theory in H[ nt : 

M$ th = {0,k'\H int 1 A t|0,fc), (63) 

Eo,k — ±1 

where |0) represents the ground FFLO state of the atoms in 
the levels |1) and |2) and the empty state of the level |3). As 
shown in Appendix [a| when uj = UJ23 + ^12/2, the transition 
amplitude is proportional to the Fourier transform p s (q x ) of 
the spin density: 

Mt, M = - 2 -^^~Ps(K-k' x ). (64) 

Then, the transition rate W$ k for the elastic scattering pro- 
cess is given by the Fermi golden rule 

Wi /M {q x ) = 2^\Mi) M \H{u k -u kl ) 
= S^ 1 ^ 1 |T 2 fc/| p 2 s (q x ) S(uj k - uo k> ), (65) 

1^12 r 

where q x = k x — k' x is the transferred momentum. The Fourier 
transform p s {q x ) of the spin density can be straightforwardly 
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FIG. 9: (Color online) Transition rate W$ k for the elastic Bragg 
scattering vs. the momentum transfer q x — k' x — k x , calculated from 
Eq. ( [66] ) for L = 121. The scattering rate has peaks at q x = 2mQ. 
The effective magnetic fields h = 1.005/i c (solid line) and h = 
1.03/i c (dashed line) correspond to the spin polarizations P « 4.1% 
and P « 4.2% in Eq. (TJ. Here we used A /£7 F = 0.2. 



computed from the Fourier expansion d58 



p s (q x ) = J iqxX p s (x)dx (66) 

oo 

= 2^ ^ p s (m) - 2mQ) + + 2mQ)] . (67) 



m=0 



We observe that the photon scattering rate Wk> ,k peaks when 
the momentum transfer q x — k x — k' x is an even integer mul- 
tiple of the momentum Q of the soliton lattice. This effect 
represents diffraction of light on the periodic modulation of 
spin density in the soliton-lattice state. An experimental ob- 
servation of this effect would give a strong evidence in favor 
of the spatially-inhomogeneous quasi- ID FFLO state. 

Equation ( [67] ) is obtained for an infinite tube with L = oo, 
but the experimental situation lfT2l corresponds to a finite sys- 
tem with L/l ~ 10 -T- 100. To check for finite-size effects, we 
calculated the integral ( [66] ) for L/l = 12 and plotted the re- 
sults in Fig. [9] for two values of h. The peaks in the scattering 
rate at q x — 2mQ are well-defined. For the experimentally 
relevant values of Q, the probe light should have the frequency 
uj in the visible spectrum in order to deliver the required mo- 
mentum transfer q x = 2mQ. This can be achieved, for ex- 
ample, by using the 2Si/ 2 —> 2P 3 / 2 transitions in 6 Li, corre- 
sponding to the wavelength of light A = 671 nm (62). Thus, 
although there may be additional complications due to con- 
finement potential, our results show that the diffraction peaks 
are observable under realistic experimental conditions and can 
give a signature of the soliton-lattice formation in quasi- ID 
fermionic superfluids. 



B. Inelastic Bragg scattering 

Now we generalize the previous discussion to the inelastic 
Bragg scattering, where the final state of the atomic system is 
an excited state \rj). The transition rate is given by the Fermi 
golden rule: 

W$ ih = 2nJ2\Mll k ( V )\ 2 8(u; kl -u k + E v - E ), (68) 



where E v — Eo is the excitation energy of the atomic system. 
The matrix element for the inelastic process is given by the 
second order-perturbation theory: 



E { 



0,k 



H 



ff in t|0,fc). (69) 



Equation 
from Eq. 



for the elastic Bragg scattering is recovered 
when uok = ojw. After some manipulations 



described in Appendix [A] we find 

4|T fc | 2 |T fc ,| 2 



S(uj k -Lj k >,k x - k' x ), 



(70) 



where 

(71) 

is the dynamical spin structure factor. Here Q = uj^ — 
and q = k x — k' x represent the energy and momentum trans- 
fers from the photon to the atoms. When the atomic system 
is in the ground state, the inelastic photon scattering is only 
possible for ft > 0, i.e., when the photon loses some energy. 
Equation ( |70| ) applies when the incoming photon energy is 
tuned to uo^ = uj 2 z + ^vi/% where the transition rate is max- 
imally sensitive to the spin-spin correlations function. For a 
different frequency, there are also contributions to the transi- 
tion rate from the density-density correlation function l62l . 

To calculate the dynamical structure factor, we use the 
imaginary-time formalism and define the spin susceptibility 



X(in n ,q)= / drdx 1 dx 2 e iQ - T - iqix2 - Xl) 



(72) 



x (T t {s z (x 2 ,t) S Z (X!,0)}), 



where Q n = 2irnT is the bosonic Matsubara frequency with 
the temperature T and an integer n. (We will take the limit 
T — >> in the final results.) In Eq. ( [72] ), we have intro- 
duced the notation s z (x,r) = ex.-p(Hr)s z (x) exp(— Hr). 
The dynamical spin structure factor S(Cl,q) is obtained from 
x(ifi nj q) by the analytical continuation: 



1 



5(0, q) — Imx(^O n — > Vt + iS : q). 



(73) 



The dynamical spin structure factor ( fTT) can be calculated ei- 
ther by direct evaluation of the averages using the exact ex- 
pressions for the Bogoliubov amplitudes similarly to the cal- 
culation in Eq. ( [5T] ) or, alternatively, by employing Green's 
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functions. Since the latter is more general and is routinely 
used in computation of various response functions, we use it 
below. The resulting expressions, of course, are independent 
of the selected method of calculation. 

The spin susceptibility x(if2 n ,g), Eq. (72), can be ex- 



pressed in terms of the normal and anomalous Green's func- 
tions for the fermions. Details of the calculations are pre- 
sented in Appendix|B] where we show that the spin correlation 
function x(i^ n ; x, x') consists of two contributions from the 
type-I and type-II processes defined below. Performing the 
analytical continuation ( [73] ), we find the spin structure factor 
in the real space: 

S{x u x 2 \ ^) = Y1 L$(x u x 2 )e-J(Sl -E^ + E v ) (74) 

+ J2^l ) (x 1 ,x 2 )e;j(n -e v - ej. 



The sums in Eq. ( [74] ) are taken over the positive energy eigen- 
values > 0, and the occupation factors 6^ are 



e 



\1V 

+ 



n F (E v - ft) - n F (E^ - ft), 



6+, = 1 - n F (E u + ft) - Tip (E^i - ft). 



(75) 
(76) 



The delta functions in Eq. ( |74| ) indicate contributions to 
S{x\,X2', fi) from two different excitation types I and II. The 
type I processes correspond to annihilation of a quasiparticle 
with the energy E v < £2 in the midgap band and creation of 
a quasiparticle in the excited state with the energy > £3 
in the upper band in Fig. [6] These processes involves spin- 
majority quasiparticles, which occupy the midgap band, so 
the factor Q~ u = 1 in Eq. ( [75] ), because n F (E v — ft) = 1 and 
n F (E^ — ft) = 0. Similar transitions for the minority spins 
are not possible, because they do not occupy the midgap band. 
The type II processes correspond to creation of two quasi- 
particles with opposite spins and the energies E v and E^. 
The minority quasiparticle can be created either in the upper 
band with E v > £3 (process Ha) or in the midgap band with 
E v < £2 (process lib), whereas the majority quasiparticle can 
be only created in the upper band with > £3, because the 
midgap band is occupied by the majority spins. In both cases, 
= linEq. (|76|), because n F (E u — ft) = n F (E v J r h) = 0. 

The matrix elements (xi, X2) and (x\, X2) for the 



type I and II processes have the typical BCS structure 1691 : 



L®(x 1 ,x 2 ) = P S^(^i) p ^^(^2) 



,b,b' 



(77) 



a, a' ,b,b' 



where the coherence factors P^Ja' ( x ) anc * TxJ^Ja' ( x ) are 
given by Eqs. (|BTT} and $BV2) . 

To calculate the matrix elements, we expand the Bogoli- 
ubov amplitudes in the Fourier series as follows 



E 



m= — 00 



u\,b(m) 
v\,b(m) 



a i(p\b+PF+Qm)x 



(78) 



Equation ( [78] ) is written for a = 1, whereas an equation for 
a = 2 can be obtained by interchanging u and v and replac- 
ing p F by — p F , as follows from Eq. p4| ). The amplitudes 
u\,b(m) and v\})(m) can be derived straightforwardly from 
Eqs. ( [27] ), ( [43] ), an d ([44] ). It also follows from the particle- 
hole symmetry, Eq. (|16|), that \u\-b{m)\ = \v\b(—m)\. Plots 
of the absolute values of several Fourier components u\^(m) 
and v\^{m) are shown in Appendix [c] 

In contrast to homogeneous superfluids, the soliton lat- 
tice breaks translational symmetry. Thus, S{xi,X2\&) and 
L^ J \xi, X2) depend not only on the relative coordinate 
x\ — X2, but also on the center-of-mass coordinate (xi+X2)/2. 
Thus, the Fourier transforms L^ 1 ^ (q, K) depend on the two 
momenta q and K corresponding to the relative and the center- 
of-mass coordinates. Substituting Eq. ([78]) into Eqs. ( |B11| ), 
( |B12| ), and ( [77] ), we obtain Eq. (B17) for the Fourier trans- 
forms L^ l \q, K). However, to obtain the inelastic scatter- 
ing rate in Eq. ( |70| ), we only need the dynamical spin structure 
factor S(Q, g) at if = in Eq. ([71]). So, we set K = in the 
rest of the calculations. 

Because of complexity of the final equations, the frequency 
and momentum dependence of S(fi, q) has to be analyzed nu- 
merically. A technical discussion is given in Appendix [C] 
Here we present a qualitative analysis and identify the under- 
lying microscopic processes in the case of a small momentum 
transfer \q\ <^p F . It follows from Eqs. (74) that 



s(n, q ) = s I (n,q) + s u (n,q), 



(79) 



where the functions Si and Sn are obtained from Eqs. ( |B17| ), 
gTg) , and ( |BT9l ) by setting K = 

Si(Sl,q)= (80) 
«i 6 > h 'i W) + ^ - E,) 

\,[jL,a,,b,b' ,{rrij } 

x ®\ii$\Pn,b' ~ Q~P\,b ~ Q{mi - mi)] <5 mi _ m ' i>m2 _ m ^, 



S u (n,q) = 

1CfJ(a,aAb',{m J })5(£l-E x -E l 

\,H,a,b,b' ,{rrij } 

x ®\p S \Pn,b' - q + Pxfi + Qim.! +m' 1 )]5 mi+m > ii 



(81) 



772-2 +m' 2 



Here the sums are taken over the energy level labels A and /i, 
the energy branch labels b and b', the Fourier indices {rrij} = 
{mi,m / 1 ;m2, m 2 }, and the Fermi points label a. The label 
a' = a in Eq. ( |B18| ) and a' 7^ a in Eq. ( |B19| ) is selected so 
that |g| <C p F . 

The numerically-calculated dynamical spin structure factor 
S(0, q) is plotted vs. the momentum transfer q and the energy 



transfer Vt in Fig. 10 for T = 0. Because the experimentally 
measurable inelastic scattering rate in Eq. ( |70| ) is proportional 
to S(fl, q), below we call S(Q, q) the signal for shortness. The 
signal is maximal at the lines A and B in the (g, O) space in 
Fig. [TUJ For a fixed q, the signal is strictly zero for the values 
of ft below the threshold represented by line B. For the values 
of Q above line B, the signal is non-zero, but maximal signal 
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band and the energy-transfer threshold represented by line B. 
While the conservation laws of energy and momentum ex- 



FIG. 10: (Color online) Contour plot of the dynamical spin structure 
factor S(Q,q) vs. the momentum transfer q and the energy-transfer 
Q. Line A represents the threshold for the type Ha processes creating 
two quasiparticles in the upper band in Fig. [6] The minimal energy 
transfer threshold ft = 2S3 is achieved at q = Q, as indicated by 
the dashed lines. Line B represents the threshold for the type I and 
lib processes, which involve one quasiparticle in the midgap band 
and another in the upper band in Fig. [6] The minimal energy-transfer 
threshold ft = 63 — £2 is achieved at q = 0. The spin structure factor 
S(fl,q) is maximal at the threshold lines A and B, as indicated by the 
bright colors. The plot is calculated for the effective magnetic field 
h — 1.05/ic corresponding to the dimensionless spin polarization 
P & 4.3% at Ao/E F = 0.2. 



is achieved at the threshold line itself. Similarly, an additional 
sharp increase of the signal is obtained when crossing line A. 

Line A in Fig. [10| originates from the type Ha processes, 
where two quasiparticles are created in the upper band in 
Fig. [6] with the energies > £3 and E v > £3. In this case, 
the minimal energy-transfer threshold is ft = 2<?3, as shown 

The type Ha pro- 



by the horizontal dashed lines in Fig. 10 



cesses are analogous to the Cooper-pair breaking in the BCS 
theory. However, the unusual feature of the FFLO state is that 
the minimal energy-transfer threshold is achieved at the non- 
zero momentum transfers q = ±Q, as shown by the vertical 
dashed lines in Fig. [TO] It is a consequence of the non-zero 
total momenta ±Q of the Cooper pairs in the FFLO state, as 
discussed in Sec.[nJ 

Line B in Fig. [10| originates from the type I and lib pro- 
cesses, which involve destruction or creation of one quasipar- 
ticle in the midgap band with the energy E v < £ 2 and cre- 
ation of another quasiparticle in the upper band with the en- 
ergy E^ > £3, see Fig. [6] Thus, the minimal energy-transfer 
threshold for line B is achieved at ft = £3 — £2, which is 
significantly lower than the BCS pair-breaking threshold of 
2Ao (roughly three times lower for the parameters used in 
Fig. [TO]). The presence of line B is a characteristic feature 
of the FFLO state. Naively, one might expect that the ex- 
cessive spin-majority fermions are unpaired, and, thus, pro- 
duce a gapless spectrum of energy excitations. However, the 
self-consistent solution discussed in Sec.[in]yields exactly one 
spin-up fermion per kink, resulting in the fully filled midgap 



plain the thresholds in Fig. 10 it is not yet clear why the sig- 
nal is maximal at the threshold lines. A technical discussion 
of this question is given in Appendix [C] Here we present a 
qualitative, heuristic explanation. Although we systematically 
wrote all equations in the paper using positive energy eigen- 
values > of the BdG equation ([14]), it is also possible 
to make a particle-hole transformation and use both positive 
and negative energies E^ . The spin susceptibility has a sim- 
pler form given by Eq. ( |B13| ) in this representation, which ef- 
fectively treats all excitations as the type I processes involv- 
ing both positive and negative energy branches in Fig. [6] In 
order to obtain a qualitatively correct result, we restrict our 
consideration to the branches with positive slope in Fig. [6] be- 
cause they correctly reproduce the normal-state energy spec- 
trum in the limit where the pairing potential vanishes. En- 
ergy and momentum conservation in interaction between pho- 
tons and Bogoliubov quasiparticles requires that the equation 
E(p — q) + ft = E'(p) is satisfied. Geometrically, the left- 
hand side of this equation represents the energy dispersion 
curve E(p) of the Bogoliubov quasiparticles in an occupied 
band, which is shifted horizontally by the transferred momen- 
tum q and vertically by the transferred energy ft. The right- 
hand side of the equation represents another, unoccupied en- 
ergy band E f (p) of the Bogoliubov quasiparticles in Fig. [6] 
If the two curves do not cross, the equation is not satisfied, 
and the signal is zero. In general, the two curves cross with 
non-parallel slopes at two points. In this case, the signal is 
non-zero, but relatively weak, because the effective overlap 
volume between the two curves at the two intersection points 
is small. However, when the threshold is approached, the two 
crossing points merge before disappearing altogether. At the 
threshold, the two curves touch with parallel slopes, which 
greatly increases the effective overlap volume between the 
two curves. Thus, we expect that the integrals of the delta 
functions in Eqs. ( [80] ) and ( [8T] ) would be greatly enhanced at 



the boundary between the domains in Fig. 10 where the con- 
servation laws can and cannot be satisfied, i.e., at the threshold 
lines A and B. Line B represents the set of values (q,ft) such 
that the displaced midgap band in Fig. [6] touches the upper 
band. It is geometrically obvious that the minimal possible 
energy transfer is ft = £3 — £ 2 at q = 0. Similarly, line 
A represents the set of values (q,fl) such that the displaced 
lower band in Fig. [6] touches the upper band. In this case, the 
minimal possible energy transfer is ft — 2<?3 at q = Q. 

These heuristic arguments indicate that the signal is en- 
hanced by maximizing the joint density of states represented 
by the delta functions in Eqs. ( [80] ) and ( [8T] ). However, these 
arguments do not take into account the structure of the matrix 
elements in front of the delta- functions in Eqs. ( [80] ) and ([8T] ). 
A more rigorous consideration is presented in Appendix |C| 

We conclude this Section by emphasizing that the dynam- 
ical spin structure factor for the soliton lattice in quasi- ID 
fermionic superfluids has much richer structure than for con- 
ventional BCS superfluids, with several important qualita- 
tive differences discussed above. Thus, the proposed inelas- 
tic Bragg scattering experiment to probe the dynamical spin 
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structure factor would help to distinguish the quasi- ID FFLO 
state from the conventional BCS superfluid in cold-atoms ex- 
periments. 



C. Radio-frequency spectroscopy 



on the ground state at T = 0, so Fermi's golden rule gives 



In Sees. |IV A\ and |IV B| we studied the case where the fre- 
quency uo of the incoming photons is detuned significantly 
from the transition frequencies between the atomic energy 
levels in Fig. [2] In that case, the fermionic atoms are only 
virtually excited from the states |1) and |2) to the state |3), 
but there are no real atomic transitions, and the photons are 
re-emitted. In this Section, we study the case where the in- 
coming photons are absorbed, and the fermionic atoms make 
real transitions from the states |1) or \2) to the state |3). Ex- 
perimentally, the state |3) in this case is typically taken to be 
the hyperfine- split state of 6 Li with F = 3/2. Let us define 
the frequency detuning ujda as 



-{ 



■ ^23 - 

■ ^23, 



^12, cr 
a 



|1>, 
|2>- 



(82) 



The interaction between the photons and the atoms is de- 
scribed by the operator ( [62] ). In this Section, we study the 
case where the photons propagate perpendicularly to the di- 
rection of the ID tubes, so the photon momentum k x = is 
zero in the x direction. Thus, we drop the index k in Eq. ( [62] ) 
and replace Tj, — » T. Then, Eq. ( |62| ) can be rewritten in the 
momentum representation, where atoms in the states \a) and 
1 3) have the same momentum p, i.e., the atoms make vertical 
transitions in Fig.[2j 



Z ' T) 



(83) 



We assume that the band |3) is empty, whereas the bands |1) 
and 1 2) are populated as shown in Fig. [2] 

The rate of transition is obtained using Fermi's golden rule 
with the perturbation H int = H a + given by Eq. ( [83] ). In 
the case where the atoms are in the normal state, the transition 
rate is 



W™TM = (2L) (2p F )\T\ 2 d(uj da ). 



(84) 



Here the factor 4Lp F represents the number of atoms in the 
occupied band. The transition rate ( [84]) i s proportional to | T | 2 , 
as opposed to |T| 4 in Eqs. ( [65] ) and (|70|) for the elastic and in- 
elastic Bragg scattering, because the latter involves absorption 
and emission of a photon in the second-order perturbation the- 
ory for the Hamiltonian ( [62] ), whereas spectroscopy involves 
only absorption in the first-order perturbation. 

In the presence of superconducting pairing, the operator 
i/j a (p) in Eq. ( [83] ) should be expressed in terms of the Bo- 
goliubov operators 7 and 7^ in Eq. ( [T3] ). Before studying the 
quasi- ID FFLO case, let us first consider a simple BCS pair- 
ing without population imbalance (h = 0). In this case, only 
the operator 7I" has a non-zero matrix element when operating 



+00 



W?$ = 2L\T\* J dp\V(pW5(Z p + ^Q + At-u da ). 

— OO 

(85) 

Here V(p) is the Fourier transform of the amplitude V(x) in 
Eq. ( [T3] ), and £ p is the normal-state energy dispersion mea- 
sured from the chemical potential in Eq. fT9] ). Given that in 
the BCS theory 



\V(p)f 



(86) 



the integral ( [85] ) can be taken be changing the variable of in- 
tegration to u = £ + y^ 2 + Aq and introducing the density 
of states p{£) = \dp/d\\/2K. Thus we obtain 1(7011711 



W^Kufr) = 2Lnp,\T\ z ( ^ ) 6{uj d(T - oj ), (87) 



where the threshold frequency is uq 
density of states p 



= Pf/V 1 + 



/~/2 . 



p. 



The 



Al)/2uj d(J pis eval- 
uated at the value of £ that satisfies the ^-function in Eq. ([S7 
and pF is the density of states at the Fermi level, see Eq. ( 48 

Physical meaning of Eq. ( [87] ) can be understood as follows. 
When the superconducting gap opens, it pushes down the en- 
ergy dispersion for the lower band, as illustrated in Fig. [3] 
So, a higher frequency is necessary to excite an atom from 
the state \a) to the state |3) compared with the normal-state 
case shown in Fig. [2] Thus, the spectral weight in Eq. ( [87] ) 
is blue- shifted to higher frequencies relative to Eq. ( [84] ). The 
detuning uoda ~ Ao corresponds to excitement of the atoms 
with \p\ ~ Pf, where p* w pp- On the other hand, at the 
threshold ujd a — > ojq, the atoms are excited from the bottom 
of the band with p — > 0, where the density of states diverges 
as ex (ujda — cjo) -1 ^ 2 . Frequency dependence of the ab- 
sorption rate ( [87] ) in the BCS state is qualitatively represented 
by the dashed curve in Fig. [TT] The spectral weight is blue- 



shifted from uoda = and monotonously decreases toward 
high frequencies. 

Now we turn to calculation of the absorption rate for the 
quasi- ID FFLO state using Fermi's golden rule 

W™ = 2^V |M(77,a)| 2 ^ a -A^, CT ), (88) 

where M(r},a) = (^,01^10, 1) is the matrix element of a 
transition from the initial to the final state, and AE v ^ a is the 
energy difference between the final and initial atomic states. 
The initial state |0,1) = |0) <E) aj|0) corresponds to the 
ground state of the atoms |1) and |2) forming the FFLO su- 
perfluid and the photon state aj|0). The final state \rj, 0) cor- 
responds to an excited state of the fermionic system with an 
atom destroyed in the state \a) and created in the state |3), i.e. 
1 77) = ipa(p)\0) ^l(p)\0), and the photon is absorbed. Us- 
ing Eq. ( [T3] ), we observe that there two kinds of the final states, 
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where a Bogoliubov quasiparticle is either destroyed 7a,ct|0) 
or created 7^ _ cr |0). According to Eq. ( [17] ), the energies of 
these states are — (E\ — ah) and E\ + ah, correspondingly. 
Because all energies are counted from the chemical potential 
/icr = ah, the energy of the atom in the third state effec- 
tively is £ p — ah, where £ p is defined in Eq. fT9] ). Putting all 
the terms together, we find 



<^°(^)=2 7 r|T| 2 2i^^ £ 

a=l 6=± £\>0 

I^S(P)| 2 n F {E x - ha) <5(£ p -E x - oj d<7 ) 
Vt b l(p)\ 2 n F {-E x - ha) d(£ p +E X - co d<7 ) 



(89) 



where the amplitudes V(p) and U(p) are the Fourier trans- 



forms of the amplitudes V^(x) and U x "' a (x) in Eq. (|78j, and 



(»), 



P = PA6 + mQ - (-l) a p F , 



^ ( 3(P) 



u\,b(m) 
v\,b(m) 



v x { %) 



v\,b( m ) 



(90) 
(91) 



Equation ( [89] ) is the main result of this Subsection. Unlike 
in the BCS theory, Eq. ( [85] ), the absorption rate in the FFLO 
state, Eq. ( [89] ), has contributions from both \U\ 2 and \ V\ 2 . The 
difference originates from the presence of the midgap band in 
the FFLO state in Fig. [6] In order to get a physical insight, 
let us consider the limiting case of the FFLO state with a van- 
ishingly small pairing potential A and focus on the majority 
fermions with a =|. In this case, | U(p) | = 1 and £ p = E\ 
for p > pp in the first term in the sum ( [89] ), so this term con- 
tributes an integral over p from pp to pp + Q/2 in Fig. [6] 
Similarly, \V(p)\ = 1 and £ p = — for p < pp in the sec- 
ond term in the sum ( [89] ), so this term contributes an integral 
over p from pp — Q/2 to pp in Fig. [6] Together, these two 
terms reproduce the contribution from the midgap band to the 
normal-state absorption given by Eq. ( [84] ). On the other hand, 
for the minority fermions with a =|, there is no contribution 
from the midgap band, because it is not populated. 

Frequency dependence of the absorption rate ( [89] ) is shown 
in Fig.fll| The dashed line shows W^^° (uJda) for the mi- 
nority fermions and is labeled as W\^ ot (uJda)^ because the mi- 
nority fermions occupy the bottom band in Fig. [3] The dashed 
line is qualitatively similar to the absorption spectrum in the 
simple BCS case. The spectrum is blue- shifted to higher fre- 
quencies, because the bottom band bends downward in Fig. [3] 
due to the superconducting gap. The contribution from the mi- 
nority fermions in Fig.[6]comes from the energy branch b = — 
with E\ > h and the Fourier component v\^(0) until the 
frequency ujda reaches the value corresponding to the quasi- 
momentum at the Brillouin zone boundary p\- = —Q/2. At 
that point, the energy conservation condition enforces a re- 
sponse from the quasi-particles in the energy branch b = + 
with E\ > h Sit the zone boundary = Q/2. To satisfy 
momentum conservation, the quasi-particles acquire an extra 
momentum — Q from the soliton lattice via the Fourier com- 
ponent v\ i+ (—l). 
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FIG. 11: (Color online) Absorption spectrum W„*3 (u>d*), 
Eq. ( [89] ), vs. the frequency detuning c^w, Eq. |82| . For the minor- 
ity fermions, a =J r , the spectrum is shown by the dashed curve and 
labeled as Wbot(^w), because it originates from the bottom band 
in Fig. [3] It is blue-shifted similarly to the absorption spectrum ( [87] ) 
for the BCS state. For the majority fermions, a =t> the spectrum 
consists of two terms Whot{^da) and Wmi&(ud<r)- The latter term 
is shown by the solid curve and originates from the middle band in 
Fig. [5] with \E\\ < hin Fig. [6] The width of this peak is of the order 
of A2, and some spectral density is red-shifted to lower frequencies. 
The plots are obtained for h = 1.25 h c , which corresponds to the di- 
mensionless spin polarization P « 5.1% assuming Ao/ Ef — 0.2. 



In contrast, the absorption rate W^\° (ud a ) for the 
majority fermions consists of two terms Whot(wda) and 
W^mid (^dcr ) • The first term comes from the bottom band occu- 
pied by the majority fermions in Fig. [3] This term is the same 
as for the minority fermions and is shown by the dashed line 
The second term VFmid(^w) comes from the mid- 
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in Fig. 

die band in Fig. [3] and Fig. [6] which is occupied only by the 
majority fermions. Because this band bends both downward 
and upward due to the opening of the gaps above and below 
it, the absorption spectrum spreads to both higher and lower 
frequencies. The term W m [d(uda) is shown by the solid line 
in Fig. 
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It is important that some spectral weight for the 
majority fermions exhibits red shifting, in contrast to the mi- 
nority fermions and the simple BCS state, which exhibit only 
blue shifting. This characteristic feature can be utilized for ex- 
perimental identification of the quasi- ID FFLO state. Similar 
theoretical conclusions were obtained in Ref. l33l . 

Absorption spectrum for the fermionic atoms was studied 
experimentally in Refs. (63] [64). The spectrum in the normal 
state consists of a single broadened peak, which is centered at 
the atomic transition frequency, i.e., at Ud a = 0. When the 
Fermi gas is cooled to degeneracy, an extra peak emerges at 
higher frequencies, followed by disappearance of the normal 
peak at the lowest temperatures [ 63 , 64 ]. Absorption spectrum 
has a characteristic threshold at the frequency ~ Aq/2E f . 
These features were interpreted as the onset of the pairing 
gap in the spectrum of single-particle excitations l63][7TH73l . 
However, we are not aware of experiments studying optical 
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absorption spectra for the quasi- ID Fermi gases with popu- 
lation imbalance in the pairing regime corresponding to the 
FFLO state, which would be very interesting. 

V. CONCLUSIONS 

In this paper, we study the quasi- ID superfluid fermionic 
condensate with population imbalance. We analyze physi- 
cal properties of this phase using the exact mean-field solu- 
tion corresponding to the soliton lattice. This mean-field ap- 
proach is valid when the inter-tube tunneling amplitude t± 
is sufficiently large, so that E F > t± > T c > t\/E F . 
We believe these conditions can be satisfied experimentally 
by choosing an appropriate optical lattice depth in the trans- 
verse direction, as well as the strength of the fermion interac- 
tion. Using the exact results for the soliton lattice, we propose 
and analyze several experiments aimed at detection of the ex- 
otic FFLO state in the quasi- ID cold-atom settings. First, we 
propose to use the optical elastic Bragg scattering to measure 
the spin-density modulation accompanying formation of the 
soliton lattice. Second, the optical inelastic Bragg scattering 
can probe the frequency and momentum dependence of the 
spin- spin correlation function and provide information about 
the quasiparticle spectrum. For both experiments, we identify 
qualitative signatures of inhomogeneous superfluidity in the 
quasi- ID FFLO state. Third, we study the difference (red vs. 
blue shift) in the absorption spectra for the radio-frequency 
spectroscopy of the majority and minority atoms. This differ- 
ence is a characteristic feature of the FFLO state in contrast to 
the conventional BCS state. Our predictions of various phys- 
ical properties of the FFLO state should help to identify this 



exotic phase in the ongoing experiments in cold atomic gases. 

In this paper, we treated the tubes in Fig. [T^ as being uni- 
form in the x direction. However, a more realistic consider- 
ation should include the effect of a confining potential along 
the tube direction, which makes the problem much more com- 
plicated. Phase diagram in the quasi- ID geometry and phase 
separation due to the confining potential were studied numer- 
ically in Ref. fljSl using a mean-field theory. The results of 
Ref. 1 16] for a fixed small t± indicate that the FFLO super- 
fluid is located at the center of the trap, similarly to the strictly 
ID case lfT2l . Therefore, we argue that the quasi- ID regime 
is the most promising for observation of the FFLO physics, 
because the atomic motion is largely one-dimensional while 
maintaining advantage of the quasi- ID Fermi surface nesting. 
At the same time, phase fluctuations of the pairing potential 
are suppressed due to the presence of weak tunneling between 
the tubes, which stabilizes the true long-range order. Phase 
diagram as a function of the interchain hopping, as well as the 
dependence of the spin polarization on t± are open questions 
at the moment and require a more systematic study. 
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Appendix A: Bragg scattering rate of light and spin structure factor for the atoms 

In this Appendix, we show how the transition rate for the Bragg scattering of light is related to the spin structure factor for 
the atoms, both for elastic and inelastic scattering. Using Eqs. ( [62] ) and ( [69] ), we obtain the following expression for the matrix 
element of transitions from the state |0, fe) representing the ground state of the atoms and a photon with the momentum k to the 
state \rj,k') representing an excited state of the atoms and a photon with the momentum k'\ 

M™ k ( V ) = J2(v,k'\Ht 1 . H„\0,k). (Al) 

The atomic level |3) is empty for both initial and final states. Using Eqs. ([62]), we rewrite Eq. in terms of the intermediate 
states \p), where the photon is absorbed, and an atom is excited from the level a to the level |3) with the momentum p 

Mj$ 9h (ri) = E / dxdJr k ,r h e- ik -*' +ik * x ( V \^)M^)\p) w 1 w (p\^ 3 (x)Mx)\0) . (A2) 

To estimate the energy difference in the denominator of Eq. ( |A2| ), we use the energy E% + p 2 /2m for the atom promoted to the 
intermediate state from the initial state with the energy E a + e(p — k x ), where we took into account momentum conservation. 
However, given that the energy corrections to the quasi-particle dispersion in the FFLO state and the momentum transfer from 
the photon are small compared with the atomic energy difference U12, we can approximately use the normal-state dispersion 
law E a + p 2 /2m for the annihilated atomic state. Then, the energy denominator becomes + E a — E$, where p 2 /2m has 
canceled out, and uj^ is the photon frequency. When the photon frequency is tuned halfway between the energy distance from 
1 3) to 1 2) and |1) (see Fig.[2j), i.e., ujk = ^23 + ^12/2, then the denominator in Eq. ( |A2| ) is equal to T^i2/2 for a =t, i- Using 
the plane-wave expansion ^s(x) = V? e%px $(p)/V2L and taking into account that (p\ffi(p)\0) = (0\^(p)\p) = 1, we find 
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from Eq. ( |A2| ) 

ny y r 

M$, k {v) = —zr-*- £ / dxdx ' e-<*' +ik * x e^ x - x 'H V \i,\(x')^(x) - ^(x')^(x)|0). (A3) 
The sum over p in Eq. ( |A3| ) gives the delta function 5(:r — which is then eliminated by integration over x' . Thus we find 

M&, h (v) = — / dxe^-^ivMxM, (A4) 

where s z (x) = \j)l(x)\j)^(x) — ^(x)^i(x) is the spin-density operator. If |^) is taken to be the ground state |0) in Eq. ( |A4| ), 
then (0\s z (x)\0) = p s (x) as in Eq. ( [50] ), and Eq. ( |A4| ) reproduces Eq. ( [64] ) for the elastic Bragg scattering. 

Substituting Eq. \AA\ into Eq. ( [68] ) and denoting q = k x — k' x ,we find the transition rate for the inelastic Bragg scattering 

^ fc = 27r 4|Tfc|2 j Tfc/|2 V / dx 1 dx 2 (0|^(x 2 )|7 ? )(7 ? |^(x 1 )|0)e^ 1 -^^(^-A^)^ (A5) 

where ^ = u k - u k > and = E v - E . Using the identity for the delta function 5(0, - AE^) = J e^-AE^t ^/2tt, the 
transition rate ( |A5| ) can be written as 

W^? fc/ =2^ 4|Tfe|2 j Tfe/|2 f / dxxdxae^c^ 1 -^) V (0|e^^(x 2 ) e -^^) (^^(^jlO) (A6) 

^12 ^-oo 27r y ~ 

= 4|Tfe|2 j Tfc/|2 / dtd Xl dx 2 e mt+i ^-^ <0|S*(*2,*)Ss(zi,0)|0} = 4|Tfe|2 J Tfe/|2 5(o; fc -a; fc ^^-^). 
cj 12 J cj 12 

Equation ( |A6| ) reproduces Eqs. ( [70] ) and ( f7T] ). 



Appendix B: Derivation of the spin-spin correlation function 

In this Appendix, we present a detailed derivation of the expression the spin-spin correlation function in Eq. ( [72] ): 

X0,r) = (f T {s z (x,r)s z (0,0)}), s z (x,t) = $\(x,t)^(x,t) - $\(x, r)^(x, r), 



(Bl) 



where Vv(^, are the fermionic field operators in the Heisenberg representation. By employing standard methods of the many- 
body theory (74), one can express the function x( x -> T ) m terms of the normal and anomalous Green's functions for the fermions, 
which we define as follows: 



g t (x,T-x',r') = -{f r {^{x,T)^\{x',r')}), g^x^x^r') = - (f T {^(a;,r)^( a ;^r , )}), 



F n (x,T-x f y) = (T t {^(x,t)^(x',t')}), 7j t (x,r;xV) = {T T {^\(x, t)^(x', t')}). 



(B2) 
(B3) 



Using Wick's theorem f74l and going into the Matsubara frequency representation, we find that the spin-spin correlation function 
dB lb can be written as a sum of two terms 



(B4) 



corresponding to the majority and minority fermions. The functions (#i, x 2 ; iCl n ) in Eq. (B4) are given by 



X 



(t) 



(xi, x 2 ; iVt n ) =T^2 [^t( x i> x ^ voo + ift n )g^(x 2 , x\\ iw) + x 2 ; itu + ^ n )J r ^(x 2 , xi; iw) 



(B5) 



X a) (^i^2;^n) = [^(^1^2;^ + in n )g i (x 2 ,x 1 ;im) +J 7 ^(x 1 ,X2iiw + iton)Ftf(x2,x 1 iizu) , 

where the sums are taken over the fermionic Matsubara frequency w. To evaluate the Matsubara sums, it is convenient to express 
Green's functions in terms of the Bogoliubov amplitudes given by Eqs. ([T3]). After some algebra, we obtain the following 
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expressions for the normal Green's functions: 



G t (x u x 2 ]izu) = 

/j,,a>,b 



inj — Efj, + h 



izu + + h 



(B6) 



izu — E jL — ft 



izu + E^ — h 



Here the summation is performed over all eigenenergies E^, the energy branch label b, and the right and left Fermi points label 
a. Similarly, for the anomalous Green's functions, we obtain 



U^jxiW^) V^i)U^{x 2 ) 
izu + En + ft 



J r n (x 1 ,x 2 ;izu) 
J r ^(xi,x 2 ;izu) 
^(x^x^izu) 
^(x^x^izu) 



izu — E^ + ft 



U$(X1)V^(X2) V^l)U^{x 2 ) 

izu + E u — h 



zzu — E fl — ft 



izu -\- E u — h 



rr( & )/ 



(B7) 



zzu — 22^ — ft 



izu — E^ + ft, izu + E^ + h 



In contrast to the conventional BCS theory, the anomalous Green's functions T^\.{x\,x 2 \izu) and ^^(xi, £2; im) are not equal 
and can be related by replacing ft — ^ — ft. Substituting Eq. ( |B7| ) into Eq. ([5]), we recover the self-consistency condition Eq. ( [26] ) 

A(x) = -gT^^{x,x]iwy^ = -2g^v x (x)u x (x) [n F (E x + ft) - n F (ft - £ A )] • (B8) 



We could have employed the particle-hole symmetry relations (16) and reduced the sums to the positive eigenenergies only. 
However, since the correlation functions (B6 ) and (B7 ) acquire the extra prefactor of two under the transformation ( 16 ), we will 
do it at the end of the calculation. 

We now proceed with the calculation of the correlation function x^\ x i, x 2\ i^n)- Substituting Eqs. (B6) and (B7) into 
Eq. (|B5|) and using the Poisson summation formula 



r E 



rip (a) — n F (b) 



n= — 00 



[mT(2n + 1) - a] [mT{2n + 1) - b] a-b 
where rip (a) is the Fermi distribution function, we obtain the following expression for \^ (#1 ? x 2 ; ifl n ): 



(B9) 



At) 



(b') 



rib'), 



n F {E IJ , - h) - n F (E„ - h) 



E V$(x 2 )V ( JJ(x 2 ) 

fi\;bb' ;aa' 



V^W^ix!) + U^ a (x 1 )U^)(x 1 ) 



- E C(*2)nJ(* 2 ) [Ul b 2(xi)VZ)(x 1 ) -V^(x 1 )Ul b J(x 1 ) 

/xA;66 7 ;aa' 



E V$( X2 )U?J(x 2 ) 



V < ^(x 1 )ul b 2(x 1 ) - U ib) a ( Xl )vli]( Xl ) 



—iQ n + — E v 




n F (-E^ - ft) - n F (-E v 


-h) 


—iVt n — E^ + E v 




n F (E^ - ft) - n F (-E u - 


h) 


-iQ n + E^ + E V 




n F (-E^ - ft) - n F [E v - 


h) 


—iVt n — E fJ — E v 



(BIO) 

The correlation function #2; ^n) can be obtained by replacing ft — > —ft. To simplify the presentation, it is convenient 

to introduce the following notation for the matrix elements: 



(Bll) 
(B12) 
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Under the particle-hole symmetry (£7 M , E v ) ->■ -£?„), one can show that («) ->■ and T^*2,(ar) -»• 

^aU*' Also ' for ^) ( -E'jU) and (E^, E v ) -> -£„), the matrix elements transform as P^^, (x) -> 

-T^l(^) and PjglUx) T( h Q ^,(x), respectively. 

Using the symmetry properties of the matrix elements, we can simplify the expressions for the correlation function. Exchang- 
ing the indices {/i, a, b} {z/, a 7 , 6'} in the expression for \^ (#1 ? #2 5 ^n) and adding it to (#1 ? ^2; *O n ), we finds 



fi A; 66'; a a ' 



^ n F (E^ - ft) - n F (E v - ft) 
—ifl n + E^ — E v 



(B13) 



Here the factor of 4 comes from combining the terms related by the particle-hole symmetry. Also, note that the sums in Eq. ( |B13| ) 
are taken over both positive and negative energies E\. Performing the analytical continuation iVt n — » ft + iS and taking the 
imaginary part of the spin- spin susceptibility as in Eq. ( [73] ), we obtain the spin structure factor in the real space 



S( Xl ,x 2 , 0) = 4 P I^IU^)P%% M [n F (E v -h)- np (E^ - h)} S(Q -E^+ E„). 

fiX;bb ; \olol' 



(B14) 



In order to understand various processes contributing the spin structure factor, it is instructive to consider the cases E^ > 0, 
E v > and > 0, E v < separately (the other cases do not contribute to S(xi , for ft, > and ^ > 0): 

5< I )(x 1 ,x 2 ,fi)=4 ]T P^,^)^ (B15) 

/2>0,^>0;66 / ;q;q; / 

S< n >(si,s 2 ,fi)=4 T^J a ,(x 1 )T^ a %,(x 2 )[l-n F (E,+ (B16) 

/i>0,^>0;66' ;aa' 

Equation ( |B15| ) describes the type I process, where a spin-majority quasiparticle is annihilated in the midgap band with the 
energy E v < ft and created in the upper band with the energy > ft. Equation ( |B16| ) describes the type II process, where two 
quasiparticles with opposite spins are created. 

We now substitute the Fourier expansions ( [78] ) into Eqs. ( |B15| ) and ( |B16| ) and perform the Fourier transform with respect to 
x\ and X2. As discussed in Sec. |IVB[ the Fourier momenta q and K correspond to the relative coordinate x\ — X2 and the 
center-of-mass coordinate (x\ +X2)/2, respectively. Thus, we obtain the Fourier transforms of the functions L^ 1 ^ appearing 
in Eqs. ([74]) and ([77]): 

Lxl fa k) = J2 , h „ x W) - - ™i + ™2 - ^2)] 

^ L — 'ai,oc' ,0,0' ,{rrij } ^ 

x 8\Pfj.,b' ~ Q~P\,b -P^p } +Pf } - Q( m i - m l + m 2 - ra' 2 )/2], 
iff & = E , fc „ , ,4? K 6, K|) - Q(mi + mi - m' 2 - m 2 )] 

x - Q+P\,b -P^p ) -P^ + + m l + m 2 + ^2)/2], 

where {m^} = {mi, m[; m/ 2 } denotes the set of indices for the Fourier components and 



(B17) 



,(b) 



=K&'), 



(B18) 
(B19) 



Using Eqs. (73), ([74), JETS) , and ( |BT9) , we obtain Eqs. ([80) and (jgT). 



Appendix C: Technical discussion of the inelastic Bragg 
scattering 

In this Appendix, we present a technical discussion of the 
elementary processes contributing to the inelastic Bragg scat- 



tering introduced in Sec. IV B 



We begin by discussing the properties of the Fourier trans- 



forms of the Bogoliubov amplitudes u\^(m) and v\^(m) de- 
fined in Eq. (78) . The absolute values of several Fourier am- 
plitudes u\^(m) and v\^(m) are plotted in Fig. 12 vs. the 
energy E\ for b = +. The particle-hole symmetry relation 
I^a, 6(^)1 = \v\,-b(— tti)\ allows one to eliminate the pro- 
cesses permitted by conservation laws but forbidden by the 
particle-hole symmetry. The plots of \v\^ (m)| vs. ft/ft c for 
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FIG. 12: (Color online) Plots of the Fourier components of the 
Bogoliubov amplitudes \u\b(m)\ and \v\b(m)\ from Eq. |78] > for 
m = 0, ±1, branch b = +, and the right Fermi point (a = 1). 
The amplitudes for the branch b = — can obtained from the ones 
shown above using \u\,b(m)\ = \vx,-b(—Tn)\. The amplitudes for 
the left Fermi point (a = 2) can be obtained using Eqs. ( [24] ). En- 
ergy dependence of these functions reflects the subtle structure of 
the wave function describing the quasi- ID FFLO state. For the pre- 
sentation purposes, we have chosen a large effective magnetic field 
h = 1.95/ic, along with k ~ 0.7 and A 2 ~ 1.4A . 



a fixed value of E\ = 1.2 £3 are shown in Fig. [13] We ob- 
serve that the Fourier components u\^(m) and v\^(m) are 
small for m ^ 0, ±1 for the experimentally relevant values of 
spin polarization, which correspond to the parameter k ~ 1/2. 
Therefore, we will only discuss matrix elements involving the 
Fourier amplitudes with m = 0, ±1 . 



Now we discuss the origin of lines A and B in Fig. 10 



for the spin structure factor S(fl,q) calculated in Sec. IV B 
First, it is useful to convert the sums over the eigenstates 
A and p in Eqs. ( [80} and ( [8T) into the energy integrations 
by introducing the density of states p(e) = J2 X S(e — E\). 
Then, one of the energy integrals can be taken by resolving 
the delta function representing the energy conservation con- 
straints En = Q ± E\ in Eqs. ( [80] ) and ( [8T] ). Another energy 
integral can be taken by resolving the delta function represent- 



ing momentum conservation: 

Si(Q,q) = Jdep b (e)pv(Sl + e)K®(e) 

x (S[p 6 /(fi + e) -q-Pb(e) + Q(m[ 
Pb(£*)Pb'(ti + £*)/C (I) (e*) 



(CI) 



mi) 



and 



S u (Sl,q) = J dep b (e)p b ^n-e)JC^\e) 

x 5\p b >(n-s) -q + p b (e) + Q(m[ +mi)] 
Pb(e*)Pb'(tt ~ e*)/C (II) (£*) 



(C2) 



where summation over a, b, b f , and {nrtj} is implied, and 
v^ 1 (e) = dp b (e)/de is the inverse group velocity of quasi- 
particles, which is related to the density of states p b (e) = 
\v^ 1 (e)\/27r. The sums in Eqs. ( |C1| ) and ( |C2| ) are taken over 
the roots e* of the equations representing momentum conser- 
vation for the type I and II processes, respectively: 



■£*) : 
- £*) 



Q- 

q 



■Pb(e*) - 
-Pb(e*) 



Q(m[ 
- Q(m[ 



mi), 
-mi). 



(C3) 
(C4) 



The roots s* exist only for the values of (q,Q) located above 
or at the threshold lines A and B in Fig. [T0[ 

The denominators in Eqs. ( |C1| ) and ( |C2| ) may vanish for 
certain combinations of the signs of b and b' when the two 
branches have equal group velocities v^ 1 = v^, 1 at the in- 
tersection point. This condition is equivalent to the condition 
discussed in Sec. |IVB| that one energy branch touches another 
one when displaced by q and ft. Vanishing of the denomina- 
tors in Eqs. ( |C1| ) or ( |C2| ) results in divergence of S(£l,q) at 
the corresponding values of (q, Q), which constitute lines A 
and B in Fig. 
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The singularity is smoothed out by a small 
but finite value of S in Eq. ( [73] ). 

Now we discuss the contributions of different processes in 
more detail. First, we consider the type I processes described 
by Eqs. ( |C1| ) and ( |C3| ), where a quasiparticle with the energy 
e < E 2 is transferred from the occupied midgap band to the 
unoccupied upper band with the energy e + > £3. Assum- 
ing that mi — m! x =0 and taking into account the energy 
spectrum shown in Fig. [6] we observe that the denominator in 
Eq. ( |C1| ) can vanish for transitions between the midgap branch 
with b = — and the upper branch with b' = +. The two energy 
branches touch for the values of (q,£l) belonging to line B in 
Fig. 10 where Eq. ( |C1| ) has singularity, and Eq. ( |C3| ) has only 
one root. The minimal energy-transfer threshold ft = S3 — 82 
is achieved at q = 0. For the values of (q, ft) located above 
line B in Fig.[l0| Eq. ( |C3| ) has two roots, and Eq. ( |C1| ) gives a 
non-singular contribution to 5(0, q). 

There are also the type I excitations with b' = b. As shown 
in Fig.[6j such processes require a momentum transfer Q from 
the lattice, i.e., m' x — mi = ±1. However, the contribution 
from these excitations is much smaller than from b ^ b' for 
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FIG. 13: (Color online) Plots of the absolute values of the Fourier 
amplitudes \v\^(m)\ vs. h/h c for E\ — 1.2 S3. The higher-order 
Fourier amplitudes decrease fast with the increase of the effective 
magnetic field h. 



the following two reasons: 1) The two group velocities have 
opposite signs, so the denominator in Eq. ([Cl]) does not vanish 
for b = b'\ 2) The matrix elements for these transitions with 



mi = ±1_ quickly decay as a function of energy e, as 

b' do not give 



shown in Fig. 12 Thus, the processes with b 
a significant contribution to S(Sl, q). 



Another contribution to line B in Fig. 10 comes from the 
type lib processes with energies e < £ 2 and ft— e > £3. Given 



that momentum conservation requires a ^ a' in Eq. (B19), 
the matrix elements for such processes involve the products 
^/x,b( m i) ^5'(mi) and Vu,b'(^i) ^,6( m i)- Similarly to the 
type I processes, the dominant contribution to the matrix ele- 
ments comes from b ^ b' and m\ = m[ =0. However, these 
processes do not give a significant contribution to S(Q, q), be- 
cause the denominator in Eq. (|C2|) does not vanish for b ^ b' 



and mi = m' x = 0. 

There is also a contribution to S(Cl, q) from the type lib 
processes with b = b' and mi — m[ — ±1. Indeed, Fig. 12 
shows that the matrix elements i/ /i+ (m / 1 = 1) u v +(m\ = Oj 
for < S2 and E\ > S3 are non-zero and strongly peaked at 
E ^ close to £2. These type lib processes with b = b' and ft = 
S3 + £2 strongly enhance 5(0, q), because the denominator 
in Eq. ( |C2| ) vanishes. Geometrically, it is a consequence of 
the peculiar nesting between the midgap and upper branches 
at q = Q and ft = £3 + £ 2 in Fig. [6] The strong enhancement 
of 5(<?3 + £2, Q) is indicated by the bright colors in Fig. 10 



Next we discuss the type Ha processes, where two quasipar- 
ticles are created in the upper band with the energies e > £3 
and ft — e > £3. These processes are possible at or above line 
A in Fig. 10 The dominant contribution comes from b' = b, 
because the denominator in Eq. ( |C2| ) can vanish in this case. 
The energy-transfer threshold is ft = 2<?3, as shown by the 
horizontal dashed lines in Fig. 10 The momentum- transfer 
threshold can be determined from Eq. ( |C4| ): 2^5(^3) = 
q — Q{m\ + m[), where we used ft = 2S3 and e = S3. 
Taking into account that | J9b (£^3 ) | = Q/2 and assuming that 
\ = 0, we find the momentum threshold at q = Q, as 



mi 



10 



shown by the vertical dashed lines in Fig. 

From momentum conservation, one might expect an exci- 
tation line starting at q = for the type Ha processes with 
b' = b, mi = 0, and m! x = —1. However, the matrix ele- 
ments , u /i ,6(m / 1 ) u^/fmi) and ^^(m^) v v $ (mi) such pro- 
cesses vanish. Figure [12] shows that the diagonal processes 
with b = b' are allowed only for mi = m' x = 0, whereas the 
matrix elements for mi =0 and rn[ = — 1 are zero. 

Finally, we consider the type Ha processes with b ^ b' . 
Given that ^+(£3) = —p-(£s) at the threshold energy, the 
momentum conservation law reads q = Q(mi+m / 1 ). Accord- 



ing to Fig. [121 the matrix elements vanish for mi 



0, 



but are non-zero for mi + m x = 1. Thus, the momentum 
threshold is q = Q for b ^ b' as well. However, the processes 
with b 7^ b' give a smaller contribution to S(Q,q), because 
the denominator in Eq. (|C2|) does not vanish in this case. 
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